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ABSTRACT 

We perform two-dimensional simulations of Alfven oscillations in magnetars, modeled as 
relativistic stars with a dipolar magnetic field. We use the anelastic approximation to general 
relativistic magnetohydrodynamics, which allows for an effective suppression of fluid modes 
and an accurate description of Alfven waves. In addition, we compute Alfven oscillation fre- 
quencies along individual magnetic field lines with a semi-analytic approach, employing a 
short-wavelength approximation. Our main findings are as follows: a) we confirm the exis- 
tence of two families of quasi-periodic oscillations (QPOs), with harmonics at integer mul- 
tiples of the fundamental frequency, as was found in the linear study of Sotani, Kokkotas 
& Stergioulas; b) the QPOs appearing near the magnetic axis are split into two groups, de- 
pending on their symmetry across the equatorial plane. The antisymmetric QPOs have only 
odd integer-multiple harmonics; c) the continuum obtained with our semi-analytic approach 
agrees remarkably well with QPOs obtained via the two-dimensional simulations, allowing 
for a clear interpretation of the QPOs as corresponding to turning points of the continuum. 
This agreement will allow for a comprehensive study of Alfven QPOs for a larger number 
of different models, without the need for time-consuming simulations. Finally, we construct 
empirical relations for the QPO frequencies and compare them to observations of known Soft 
Gamma Repeaters. We find that, under the assumptions of our model and if the magnetic field 
of magnetars is characterized by a strong dipolar component, and QPOs are produced near the 
magnetic pole, then one can place an upper limit to the mean surface strength of the magnetic 
field of about 3 - 8 x W^G. 
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1 INTRODUCTION 

The phenomenon of Soft Gamma Repeaters (SGRs) may allow us 
to determine fundamental properties of strongly magnetized, com- 
pact stars in the near future. SGRs produce giant flares with peak 
luminosities of 10^"* - 10''® erg/s, which display a decaying X- 
ray tail for several hundred seconds. Up to now, three such gi- 
ant flares have been detected, coming from SGR 0526-66 in 1979, 
SGR 1900+14 in 1998, and SGR 1806-20 in 2004. Quasi-periodic 
oscillations (QPOs) have been observed in the X-ray tail of the 
latter two sources, following the initial discovery by llsrael et al] 
l l2005h (seelWatts & Strohmaver (2006) for a recent review). The 
timing analy sis of these QPOs allowed the determination of their 
frequencies l lWatts & Strohmaveill200^ . which are approximately 
18, 26, 30, 92, 150, 625, and 1840 Hz in the SGR 1806-20 gi- 
ant flare, and 28, 53, 84, and 155 Hz in the SGR 1900+14 gi- 



* E-mail:cerda@mpa-garching.mpg.cle 
t E-mail:mksterg@auth.gr 
X E-mail:j. antomo.font@uv.es 



ant flare. The frequency of many of these oscillations is similar 
to what one would expect for torsional modes of the solid crust of 
a compact star. This observation is in support of the proposal that 
SGRs are magnetars, i.e. compact objects with very strong mag- 
netic fields (Duncan & Thompson 1992). During an SGR event, 
torsio nal oscillations in the solid crust of the star could be ex- 
cited l lDuncanlll998h . leading to the observed frequencies in the 
X-ray tail. 

While a lot of theoretical work initially foc used on the 
torsiona l oscillation s of the soli d crust (see ISotani et al.l 
( l2007ah : ISamuelss on & AnderssonI ( |2007|) and refe rences 
therein), more recently elasto-magnetic oscillations jLevinl 
l2006l : lGlampedakis et alfcood: lLevinil20a7l : lLe3l2007ll2008l) and 
pure Alfven oscillations jSotani et al .120081) have become attractive 
in the context of the observed oscillations in SGR giant flares. In 
particular, Leviri ( 200^ and [oiampedakis et al. (2006) stressed 
the i mportance of crust-core coupling by a global magnetic field, 
while !Leviiil ( l2007l) presented a toy-model of the Alfven continuum 
and demonstrated that long-lived QPOs can be pro duced at the 
edges or turning points of the continuum. In addition, ISotani et al.l 
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( l2008h . using two-dimensional linear simulations, found that in a 
magnetar model that uses ideal magneto-hydrodynamics (MHD) 
with a microphysical equation of state, long-lived QPOs appear 
near the magnetic axis and within the region of closed magnetic 
field lines. In addition, they found that an infinite number of 
"overtone QPOs" exist, with frequencies at integer multiples of the 
fundamental frequency, in rema rkable agreement w ith observed 
QPO frequencies. Furthermore, ISotani et alj ( 12008 ) found that 
all Alfven QPO frequencies are described by simple empirical 
relations, valid for any equation of state and magnetar mass, 
initiating the formulation of magn etar asterose ismology based on 
global Alfven QPOs (see also Samu elsson & A ndersson (2007) for 
a formulation of asteroseismology based on torsional modes of the 
crust). 

The Alfven QPO model is attractive, since it reproduces the 
observed near-integer ratios of 30, 92 and 150 Hz in SGR 1806- 
20. Since these QPOs appear at odd-integer ratios, our current re- 
sults suggest that they could correspond directly to the odd-parity 
QPOs appearing near the magnetic axis. In addition, the Alfven 
QPO model is attractive for the reason that it is difficult to explain 
all of the frequencies of 18, 26 and 30 Hz in the SGR 1806-20 giant 
flare as being due to crustal shear modes, since the actual spacing 
of torsional oscillations of the crust is larger than the difference be- 
tween the 26 and 30 Hz frequencies. Only one of these frequencies 
could be the fundamental, £ = 2, m = torsional frequency of the 
crust, as the first overtone has a much higher frequency. Note that 
the 18 Hz frequency could be explained ei ther by the Alfven QPO 
model or as the fundamental crust mode l lSteiner & Wat3ll2009l) . 
Similarly, the spacing between the 625Hz and a possible 720Hz 
QPO in SGR 1806-20 may be too small to be explained by con- 
secutive overtones of crustal torsional oscillations. The spectrum 
of observed QPOs may thus include both crustal normal-mode fre- 
quencies and Alfven QPO frequencies. 

Here, we present the first two-dimensional numerical sim- 
ulations of Alfven oscillations in magnetars in the anelas- 
tic approximation to g eneral relativistic magnetohydrodynamics 
dSonazzola et al.ll2007 ). in which fluid modes are suppressed, ob- 
taining efficient and accurate simulations of Alfven waves. In ad- 
dition to the numerical simulations, we also compute Alfven os- 
cillation frequencies along individual magnetic field lines with a 
semi-analytic approach, employing a short-wavelength approxima- 
tion. The continuum obtained using standing-wave oscillations in 
the semi-analytic approach agrees remarkably well with QPOs ob- 
tained via our two-dimensional simulations, allowing for a clear in- 
terpretation of the QPOs as corresponding to turning points of the 
continuum. This agreement will allow for a comprehensive study 
of Alfven QPOs for a large number of different models, without 
the need of time-cons uming simulations. In addition to the QPOs 
found in I Sotani et alj (12008) (which were antisymmetric with re- 
spect to the equatorial plane) we also find a separate sub-family of 
QPOs that are symmetric. Furthermore, our implementation of the 
zero-traction boundary condition at the surface of the star shows 
that the QPOs generated near the axis have a nonzero amplitude 
at the surface, which would allow these oscillations to reach the 
magnetosphere. 

Our semi-analytic approach is not restricted to a dipolar field, 
but can be extended to other magnetic field configurations, while 
our anelastic approximation does not suffer fr om the nume r ical in - 
stabilities encountered in the linear study of ISotani et al. I j2008t) . 
Furthermore, our anelastic approach could be used to extend the 
current sudy by including shear waves in the crust, which may 
prove to be important for the interpretation of observed QPOs. 



The paper is organized as follows: Section 2 describes the the- 
oretical framework of the general relativistic MHD equations in the 
anelastic approximation. Section 3 discussed the equilibrium mod- 
els of magnetars we consider. The Alfven continuum and the condi- 
tions for standing waves along individual field lines are studied with 
a semi-analytic approach in Section 4. In Section 5 we present our 
axisymmetric simulations in the anelastic approximation, while in 
Section 6 we present empirical relations for the obtained QPO fre- 
quencies and compare them to observations of known SGR sources. 
Finally, we discuss our present results in Section 7. Unless stated 
otherwise, we adopt units of c = G = 1, where c and G are the 
speed of light and the gravitational constant, respectively, and the 
spacetime metric signature is ( — , +, +, +). Latin (Greek) indices 
run from I to 3 (0 to 3) and boldface symbols indicate vectors. 



2 THE EQUATIONS OF GENERAL RELATIVISTIC MHD 
IN THE ANELASTIC APPROXIMATION 

In our theoretical framework we use the 3+1 split of spacetime in 
which the metric reads 



~a dt^ + jijidx' + f3' dt)(dx^ + (5^ dt), 



(1) 



where a and /3* are the lapse function and the shift vector, respec- 
tively. Furthermore, we adopt a conformally flat approximation in 
which the spatial 3-metric is 

,4 - 



(2) 



where (fi is the conformal factor and "fij denotes the flat metric. The 
energy-momentum tensor of a magnetized perfect fluid is 



r"" = {ph + b^) u^u" +[p+"— 



(3) 



where p is the rest-mass density, h is the specific enthalpy, P is the 
isotropic fluid pressure, is the 4-velocity of the fluid, f is the 
magnetic field measured by a comoving observer, and g^j^ is the 

4-dimensional metric tensor . 

Following lAnton et al.l ( |2006|) we define the following set of 
conserved quantities (which are generalizations of the rest mass 
density, momentum density and energy density, respectively) 



D = pW, 

S^ = {ph + b'')W^v,-aUb'\ 



T = (ph + b'')W'^- 



-a^ib^'f -D, 



(4) 
(5) 

(6) 



with being the 3-velocity of the fluid and W = avP the 
Lorentz factor. This choice allows to cast the general-relativistic 
magneto-hydrodynamics (GRMHD) equations in the form of a 
flux-conservative hyperbolic system of equations, namely 



dt 



(7) 



where y/^-g = 7 ~ <iet(7ij) and where the state vector, 

the flux vector, and the source vector are given, respectively, by 



U = [D,S„B% 



v'B" 



W 



(8) 



(9) 
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Table 1. Magnetar equilibrium models. From left to right the columns report the model name, the circumferential equatorial radius R, the 
coordinate equatorial radius r^, the ratio of polar to equatorial coordinate radii (minus one), the gravitational mass M, the magnetic field 
strength at the pole Spdar, the central rest mass density p^, the current density jq, and the form of the current density in the interior (for 
model SI the value of j'q was not available). Model MNS2 (in bold) is our reference model. 



Model 


ij[km] 


T'e [km] 


rp/r'e - 


1 


M [Mq] 




Pc[gcm 3] 


jo [Am 2] 


current function 


MNSl 


14.155 


11.998 


8 X 10- 


6 


1.40 


6.5 X 1014 


7.91 X 10" 


2 X 10^3 


Bocauet et al. (1995) 


MNS2 


14.157 


11.999 


8 X 10- 


4 


1.40 


6.5 X lO^^^ 


7.91 X 10" 


2 X 10" 


Bocauet et al. (1995) 


MNS3 


14.168 


12.006 


5 X 10- 


3 


1.40 


1.6 X 


7.91 X 10" 


5 X 10" 


Bocauet et al. (1995) 


LMNS2 


15.153 


13.322 


10-3 




1.20 


5.1 X 10^5 


5.47 X 10" 


2 X 10" 


Bocauet et al. (1995) 


HMNS2 


12.444 


9.941 


4 X 10- 


4 


1.60 


8.4 X 10^5 


1.37 X lOi""' 


2 X 10" 


Bocauet et al. (1995) 


SI 


14.153 


11.912 







1.40 


4.0 X lO^s 


7.91 X 10" 


N/A 


Sotani et al. (2007a) 



^ yA"^ dgij,v 



0,0,0 



(10) 



where = — P"" /a. 

In the above expressions, denotes the magnetic field com- 
ponents meas ured by an Euleri an observer . The relation between 
6^ and is jAnton et al.ll2006h 



6^ = 



W 



(11) 



Notice that we consider a barotropic equation of state (EOS), 
P = Kp^, for which an evolution equation for the energy r is not 
required. This type of EOS is in particular appropriate for studying 
those small-amplitude stellar oscillations for wh ich thermal effects 
can be ignored. As shown in lAnton et"al] ( 12006 ) the eigenvalues 
associated with the Jacobian matrices of the GRMHD system. A*, 
are an explicit function of the metric tensor, the fluid and magnetic 
field variables and the speed of sound, Ca, i.e. 



(12) 



The dependence on the speed of sound is a consequence of the 
pressure terms appearing in the fluxes which lead to the presence 
of a term dP/dU in the Jacobian matrices dF^ /dU . For the case 
of a barotrope this term reads 



(13) 



dP _ dP_dp_ _ 2dp_ 
dU ^ dpdU ^ dU ■ 

The m ain idea of the anelastic approximation dSonazzola et al.l 
l2007h is to eliminate sound waves, in order to reduce the time- 
step restrictions of numerical codes that are based on time-explicit 
schemes. In the presence of sound waves such restrictions are im- 
posed by the dependence of the eigenvalues of the hyperbolic sys- 
tem on the speed of sound. Consider the following modified equa- 
tions in which the pressure terms in the fluxes have been moved to 
the sources: 



w 



,v B 



V B 



'2 dx^ 



,0,0,0 



(14) 



(15) 



With this choice, the speed of sound does not appear in the eigen- 
values of the corresponding Jacobian matrices. In fact, these modi- 
fied eigenvalues can simply be obtained from the original eigenval- 
ues by taking the limit of vanishing speed of sound, i.e. 



\\-f^j,U)=\\-i.J,U,C.^Q). 



(16) 



With this modification the time-step required to satisfy the stability 
criterion imposed on hyperbolic equations by the Courant condi- 



tion is much less restrictive, since it is not affected by sound waves. 
We note however, that this procedure cannot be applied in a more 
general context, since the source terms may become stiff due to the 
inclusion of the pressure, which may lead to unstable numerical 
evolutions of the system. Nevertheless, in the present anelastic ap- 
proximation this approach is valid since we further assume that the 
pressure remains constant in ti me. In the case of a barotrope, this 
leads to teonazzola et al.|[2007h 



dt 

which implies 



-gDff 



dx^ 



= 0. 



(17) 



(18) 



The anelastic approximation can be used in our numerical study 
since we deal with low-amplitude torsional oscillations of an equi- 
librium star, which are of axial parity. Axial velocity perturbations 
couple only weakly to density perturbations, so that the density 
(and pressure) can be considered constant. 

In order to stabilize the numerical code, we subtract from the 
source terms of the momentum equations the initial equilibrium 
value of the flux-source balance (which should be zero but it is 
not due to discretization error). For axial oscillations at the linear 
level, only the perturbation in B"^ and 5^ are nonvanishing, so that 
we can keep B' , B®, Sr, and Se constant in time (as long as the 
oscillations have small amplitude). 



3 EQUILIBRIUM MODELS 

In our analysis we use two different methods to build rela- 
tivistic equilibrium models for non-rotating magnetars. The first 
method (valid up to moderate magnetic field strengths relevant for 
magnetars) relies on first computing a (nonmagnetized) Tolman- 
Oppenheimer-Volkoff (TOV) solution, and then separately com- 
puting the magnetic field structure. Since the TOV solution is de- 
scribed in Schwarzschild coordinates we make the transformation 
to quasi-isotropic coordinates. We further assume that the con- 
figuration of the magnetic field is that of a pure dipole. Details 
on the numerical method for constructing the magnetic field, as 
well as a representative figur e of the magnetic field lines, can be 
found in Sotani et al.' ( 12007 ah . The same configuration was used 
in Sotani et al. (2008). 

The second method we use for constructing equilibrium mag- 
netar models is based on the self-consistent solution of magnetized 
stars with a dipolar magnetic field, where all the effects of the 
magnetic field on the matter and spacetime are taken into account. 
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presented bv lBocguet et alj il995h . These equilibrium models are 
computed using the LORENE library 0. We construct non-rotating 
polytropic equilibrium models with F = 2 and K = 1.455 x 10^ 
(cgs units). The specific central enthalpy is chosen to be he = 
1.256 and the form of th e current densi t y whi ch gives rise to the 
magnetic field is given in lBocguet et al.l jl995l) . Table [T] describes 
the main properties of all equilibrium initial models used in our 
study. The equilibrium model MNS2 is used in the following Sec- 
tions as a reference model, with a mass of 1.4i\f0 and a magnetic 
field strength at the pole equal to -Bpoic ~ 6.5 x 10^^ G. 



4 A SEMI-ANALYTIC MODEL 

The time evolution of the torsional Alfven oscillations in the neu- 
tron star interior is a computationally expensive task, even in the 
anelastic approximation defined above. This fact restricts the num- 
ber of models that can be studied with our numerical approach. 
Therefore, in addition to the numerical GRMHD simulations that 
will be presented in later Sections, we have developed a simpli- 
fied approach for computing the oscillation spectrum. This semi- 
analytic model is based on the computation of Alfven wave travel 
times along individual magnetic field lines (assuming a short- 
wavelength approximation) and requires only information on the 
initial equilibrium model and magnetic field configuration. 

In the linear regime and in the limit of short wavelengths 
(compared to the length of magnetic field lines inside the star) tor- 
sional Alfven oscillations can be treated as axial perturbations trav- 
eling strictly along individual magnetic field lines. This is evident 
in the local characteristic structure of the general relativistic MHD 
equations. The Alfven eigenvalues of the system of equations for a 
given direction n is ( Anton et aL .2006) 

6 ■ n ± \/C 11 ■ n 



A 



(19) 



where C = ph + and u is the fluid speed. For a non-rotating 
initial equilibrium configuration = 0, and therefore the eigen- 
value has a maximum when n is parallel to the magnetic field b. 
Then, we can construct the Alfven velocity vector as 

ab 



(20) 



i'a(a;), 



(21) 



which describes two possible waves along a magnetic field line, 
traveling in opposite directions. 

In this approximation, the path of an Alfven wave inside the 
star can simply be found by integrating the equation 

dx 
It 

starting from a given initial location xq inside the star Since the re- 
sult of this integration coincides with a magnetic field line passing 
through xo, we refer hereafter to the path of the Alfven waves sim- 
ply as magnetic field lines. We integrate l l21t numerically using a 
second order Runge-Kutta method. Due to symmetry, only the part 
of the magnetic field lines in the upper plane is computed.The start- 
ing points for this integration are located in the equatorial plane, 
and the integration ends either at the surface of the star (for open 
magnetic field lines) or at another point in the equatorial plane (for 
close magnetic field lines). 

Since the topology of the magnetic field lines considered here 

^ http://www.lorene.obspm.fr 
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Figure 1. Schematic illustration of the integration domain. In the region 
of open magnetic field Hnes (pink area) lines start at the surface (point a), 
cross the equatorial plane (point a') and end at the opposite end of the sur- 
face (point a"). For closed magnetic field lines (blue area), points b and b" 
coincide (left panel). The same magnetic field lines and points are shown in 
coordinates (x, §), which are adapted to Alfven waves traveling along field 
lines (right panel, see text for details). 



is rather simple, we can parametrize all magnetic field lines inside 
the star by their radial position r in the equatorial plane, normalized 
to the radial position Vc of the location in the equatorial plane where 
close magnetic field lines have zero length (see Fig[Tll. We call this 
parameter x = f/Tc, ranging from to 1 (notice that field lines 
ending in the region < r < re in the equatorial plane, where 
is the equatorial radius, originate in the region r < r^). Using 
Alfven waves that travel along individual magnetic field lines as 
tracers (starting from the equatorial plane), the location of a spe- 
cific point along a magnetic field line is a function of the parameter 
X and of the arrival time ta of an Alfven wave at this location i.e. 
X — x{x;ta)- Furthermore, we indicate as ttot(x) twice the to- 
tal travel time of an Alfven wave traveling along a magnetic field 
line characterized by the parameter x> starting from the equato- 
rial plane and ending either at the surface of the star or at another 
place in the equatorial plane. Then, the dimensionless parameter 
^ = ta{r,9; x)/itot (x) ~ 1/2, characterizes the location of points 
along an individual magnetic field line x- For open magnetic field 
lines inside the star (see Fig.[T]l, ,^ takes the value at the lower 
end of an open magnetic field line (i.e. at the surface of the star, 
below the equatorial plane) and the value of 1 at the other end (at 
the surface of the star, above the equatorial plane), while for closed 
field lines ^ can take the value of at any point and the value of 1 
at the same point, after completing a full cycle. In this way, we can 
transform a magnetic field line from polar coordinates (r, 9) to the 
dimensionless magnetic-field-line-adapted coordinates (x,C)- The 
above procedure could also be extended to more complicated or 
to three-dimensional topologies by adding more parameters and/or 
regions. 

Next, assume that x\ ^) is a displacement that describes 
traveling waves along a specific magnetic field line x- Since we 
focus on individual field lines, we can simplify notation and only 
show explicitly the dependence on ^ and t in what follows. For a 
traveling wave, Y satisfies trivially the wave equation 



= a 



(22) 



where a is the wave speed. It is easy to show that, owing to our use 
of the dimensionless arrival time ^ — ta/ttot for characterizing 
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the location of points along a magnetic field line, the wave speed 
a is constant and equal to ct = 1 /ftot • In other words, the depen- 
dence of the Alfven speed on position has been absorbed into the 
definition of 5. 

To find a correspondence with QPOs, we are interested in so- 
lutions in the form of standing waves 



y(C,t) = A(C)cos(27r/t + ^), 



(23) 



where / is the oscillation frequency and <^ is a phase. The spatial 
dependence can be written as 



A(^) = a sin(K^ + I 



(24) 



where n is the wavenumber, (f> a spatial phase, and a the amplitude 
of the oscillation. Substituting l |23l l into l l22t . the resulting disper- 
sion relation is simply 



(25) 



If we assume the existence of standing waves for a single magnetic 
field line, then t) — Y{^, t + tcydc), where tcycic is the time 
it takes for an Alfven wave to complete a full cycle along the mag- 
netic field line and return to its initial position (for open magnetic 
field lines, the wave will be reflected at the surface above and be- 
low the equatorial plane before completing the cycle). An infinite 
number of discrete oscillation frequencies corresponds to standing 
waves 



fn 



cycle 



n = 1,2,3., 



(26) 



with /i = 1 cycle being the frequency of the fundamental mode. 

To simplify further the analysis, we decompose the spatial part 
of the perturbations into odd (-) and even (-I-) parity components 
with respect to the midpoint ^ = 1/2 



(27) 



If the background equilibrium (unperturbed) model is symmetric 
with respect to the equatorial plane, then ^ = 1/2 coincides with 
the equatorial plane, and the two components of the perturbation 
are just symmetric (-I-) and antisymmetric (-) oscillations with re- 
spect to the equatorial plane. Note however, that such a split would 
also exist in a general case, without equatorial plane symmetry, 
since there is always a midpoint ^ = 1/2. The symmetry of the 
oscillations fixes the spatial phase of each component, so that 



AiiO = cosK(C-l/2)], 
A'iO = an sin[«;„(e-l/2)]. 



(28) 



This decomposition can be related to Eq. ( I24l > for each frequency 
fn , considering 



+ a„ 



(j)n = arctan (a„ /a^) — «:„/2. 



(29) 



To completely determine the solution of the problem one has 
to further consider the topology of the magnetic field lines and the 
boundary conditions. Due to the solenoidal condition satisfied by 
the magnetic field, only two possible magnetic field line topologies 
are possible for any field inside the star: open or closed field lines, 
which we consider separately, next. 



Y(i,t) = y(o,t). 



(30) 



Therefore, the time that takes for a perturbation to complete a cycle 
is just the time to travel along the line, fcycic ~ ttot- The frequency 
of each standing wave is 
n 

fn = — ; n=l,2,3... (31) 

ttot 

Using the dispersion relation ( I25t the wave number for each 
standing wave can be found as 



Kn — 2tv n 



n = 1,2,3, 



(32) 



Therefore, only standing waves with discrete frequencies /„ are 
possible along a single magnetic field line. The number of nodes 
for each standing wave is just 2 n, and, using l l28t the symmetric 
and antisymmetric component of the oscillations are 



cos(27rn5), 
sin(27rn^). 



(33) 



If we consider purely symmetric oscillations, then a~ = and 
therefore the spatial phase is fixed to (j)„ — 0, for every n. In 
the case of purely antisymmetric oscillations a+ = and hence, 
4>n ~ 7r/2. In the general case, without any symmetry in the per- 
turbations, the spatial phase depends on the initial condition and 
can be different for each n. 

Since in this region of closed field lines each magnetic field 
line is decoupled from the others, the frequency forms a con- 
tinuum for each overtone. 



4.2 Open magnetic field lines 

In this case, each magnetic field line crosses the surface of the star 
at two points, ^ = and ^ = 1. The time it takes an Alfven wave 
to complete a cycle, i.e. to travel along the magnetic field line and 
return back to the starting point is tcycic ~ 2 ttot ■ Therefore, the 
allowed frequencies for standing waves are 



fn 



2ttot 
and therefore 

K„ = TT n : 



1,2,3., 



n = 1,2,3., 



(34) 



(35) 



The symmetric and antisymmetric components read in this case 



A-niO = an sin (tth^) ; n=l,3,5.. 

^n(0 = COS (7rn^) ; n — 2,4,6.. 

A-niO ~ dn cos(7rn^) ; n=l,3,5. 

= a;;sin(7rO ; n = 2,4,6.. 



(36) 
(37) 
(38) 
(39) 



It is important to note that symmetric oscillations with even n and 
antisymmetric oscillations with odd n, i.e. Eqs. J37t and l l38t , fulfill 
the conditions yln (0) — An{l) = at the boundaries. On the other 
hand symmetric oscillations with odd n and antisymmetric oscil- 
lations with even n i.e. Eqs. l l36t and l l39t satisfy the conditions 
d^An{0) = d^An{l) = at the boundaries. Therefore, the pos- 
sible standing waves depend on the symmetry of the perturbation 
with respect to the equatorial plane and on the boundary conditions 
imposed at the surface. 



4.1 Closed magnetic field lines inside the star 

We consider first magnetic field lines that close inside the star. In 
this case, standing-wave perturbations fulfill spatial periodicity 



4.3 Boundary conditions at the surface 

The boundary conditions at the surface of the star should ensure 
the continuity of traction. In the present approach, we neglect the 
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0. 



(40) 



dynamics in the magnetosphere, so that the traction is zero at the 
surface. If we do not consider the presence of a crust, the zero trac- 
tion condition for Alfven oscillations at the surface would lead to 
the boundary condition that the gradient of Y along each field line 
should vanish at the surface. However, in reality there exist a solid 
crust in the star, which occupies a region of about 5 to 10% of the 
radius. This changes the boundary condition at the surface of the 
star. Furthermore, the oscillations behave differently in the fluid re- 
gion, where only Alfven waves are possible, and in the crust, where 
both Alfven and shear waves appear. In the case that the magnetic 
field is not too strong, the Alfven speed in the crust is comparable 
or smaller to the shear speed. In this case we can make the ap- 
proximation that the boundary condition at the surface of the star 
corresponds to the zero-traction condition for shear waves in a non- 
magnetized star, i.e. 

dY 
dr 

Furthermore, we neglect the shear waves in the crust. This approx- 
imation is valid in two cases. First, if shear and Alfven waves have 
similar velocities, then the velocity of the perturbations in the crust 
is underestimated (but correct in order of magnitude), and the di- 
rection of propagation is not exactly along the magnetic field line. 
Second, if shear waves dominate the dispersion relation at the crust 
(low magnetic field), then the perturbations traveling along mag- 
netic field lines inside the core and reaching the crust interface, 
will propagate almost immediately inside the crust (in comparison 
with the travel time inside the core) to reach the surface where the 
boundary condition assumed in this work holds. Therefore our ap- 
proximation in this limit has the meaning of considering a very thin 
crust, and hence neglecting its inertia. In both of the cases, the ap- 
proximation can produce small quantitative differences in the fre- 
quency spectrum and damping time, but the qualitative behavior 
should not change. This approximation should be valid as long as 
the magnetic field is not too strong (possibly larger than 10^'' G). 
This boundary condition, contrary to the case in which the presence 
of a crust is completely neglected, introduces a coupling between 
the different magnetic field lines which affects the time evolution 
of the QPOs that appear at turning points. 

Since the perturbation is a function xi ^) (in the follow- 
ing we include again the explicit dependence on x) the boundary 
condition l|40l> becomes 



(41) 



200 



dr dx dr 
Using the standing-wave solution ( 123b yields 



cos(27r/ 1 



-Ai^,x)sm{2nft + ^) [ 2iit-^ + 



dA{i,x) 
dx 

df 



cos(27r/ 1 + (fi) 
dx 



dip 
dx 



dr 



0. (42) 



Due to these boundary conditions at the surface, standing-wave so- 
lutions are possible in two different cases. For perturbations which 
vanish at the surface, i.e. for solutions of the form l l37t and l l38t . the 
only solution of ( 142b is the trivial solution for vanishing amplitude. 
On the other hand, for solutions of the form ( |36b and J39b , (sat- 
isfying d(_An = at the surface) a standing-wave solution exists 
if 



dx dx) dr 



da^dx ^ Q 
dx dr 



(43) 
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Figure 2. Standing-wave continuum and QPOs for the reference model 
MNS2. Filled circles denote the (upper and lower) QPOs. See text for de- 
tails. 



dl 

dx 



dip 
dx 



da 
dx 



(44) 



i.e. if f{x), fix) 'tnd a(x) have a common turning point, or if they 
are constant. The second case is when 



^=0, 
dr ' 



(45) 



i.e. for those magnetic field lines that cross the stellar surface per- 
pendicular to it. In axisymmetry, this always holds for the mag- 
netic field line along the pole. Therefore, only at the magnetic pole 
an exact standing wave solution is possible. The spatial part of the 
solution is 



AniO = cos(-7rn5), 



(46) 



which has n nodes inside the star. As we discuss next, one can still 
define a standing wave along any magnetic field line, on a short 
enough timescale. 



4.4 Standing-wave continuum and QPOs 

According to the standing-wave solution for the region of open field 
lines, presented in the previous section, as long as 



2 A « 1, ^ 
dx dx 



0, 



|^~0. 

dx 



(47) 



the time-evolution of an initial perturbation will still be close to 
a standing wave. Therefore, the characteristic timescale on which 
phase mixing occurs and a standing wave is damped, is 



Td 



2^ 



(48) 



These conditions can be fulfilled only in two cases. First, if 



Consequently, on a timescale smaller than Td, one can define a 
standing-wave continuum for magnetic field lines throughout the 
star, for each overtone. 

In Fig. |2] we display the different branches of the standing- 
wave continuum of frequencies for the reference model MNS2, 
in a frequency vs. magnetic field line plot. The continuum can 
be separated into two parts, one for the region of open field lines 
(0 < X < 0.66) and another for the region of closed field lines 
(0.66 < X < l-O)- Standing waves in the region of open field lines 
can be divided into two sub-families, according to their symmetry 
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Table 2. QPO frequencies for the reference model MNS2, in the 
anelastic approximation. In parenthesis we give the relative differ- 
ence with respect to the semi-analytic approach. The frequency res- 
olution of the numerical simulations is 1.25 Hz. 



"upper" QPOs f (Hz) "lower QPOs" f (Hz) 





35.1 (4.%) 




41.30(1.0%) 




72.6 (0.05%) 


Li 


83.85 (0.6%) 




110.2(1.1%) 


Li 


126.4(1.0%) 




147.8(1.8%,) 


Li 


170.2 (2.0%) 




182.8(0.7%)) 


Li 


210.3 (0.9%) 


^2 


220.4(1.2%,) 


Li 


251.6(0.6%) 




256.7(1.0%) 


Li 


295.4(1.2%) 


^3 


294.3(1. 3%c) 


Li 


336.7 (0.9%) 



with respect to the equatorial plane, with distinct frequencies. In the 
region of closed field lines, these two sub-families have identical 
frequencies. Notice that, for each overtone, the branch which cor- 
responds to symmetric standing waves in the region of open field 
lines smoothly joins to the two overlapping branches of symmetric 
and antisymmetric standing waves in the region of closed magnetic 
field lines. 

For the particular magnetic field configuration used here, the 
standing-wave frequencies have a local maximum at x = (i.e. 
at the pole). Notice that conditions ( |47| ( are always fulfilled at the 
pole, since it is a turning point of the continuum (due to axisymme- 
try). Hence, near the pole, we expect local standing waves to give 
rise to long-lived QPOs. The predicted frequencies of these QPOs 
at the pole (the "upper" QPOs) are shown in Fig. ^ and labeled as 
and Ui~\ for symmetric and antisymmetric standing waves, 
respectively. Here, the index n labels the different overtones (start- 
ing with the fundamental frequency, n = 0) and also coincides 
with the number of nodes appearing between center and pole. 

In the region of closed field lines, the standing-wave contin- 
uum for each overtone features a local minimum. Notice that the 
continuum for the most part of this region is nearly flaQ implying 
a slow rate of phase-mixing. These local minima also give rise to 
long-lived QPOs in our numerical simulations (the "lower" QPOs), 
which are labeled as L^i^ (the frequency of symmetric and anti- 
symmetric standing waves coincide). Here, the index n again labels 
the different overtones (starting with the fundamental frequency, 
n — 0). For the fundamental frequency, two nodes appear along 
a closed field line, corresponding to a single nodal line crossing 
the region of closed field lines. For each consecutive overtone, our 
semi-analytic model predicts an additional nodal line. 



5 NUMERICAL SIMULATIONS 

Our second approach in studying Alfven QPOs in magnetars in- 
volves two-dimensional num erical simulations u s ing th e nonlinear 
GRMHD code described in ICerda-Duran et alj ( l2008h . in which 
we implemented the anelastic approximation presented in Sec. |2] 
We use spherical polar coordinates under the assumption of ax- 
isymmetry. This code solves the GRMHD equations coupled to 

^ At X = 1 the frequency tends to infinity, because the length of the mag- 
netic field line tends to zero. 



the Einstein equations for the evolution of a dynamical spacetime 
(in the approximation of spatial conformal flatness) and has been 
developed with the main objective of studying astrophysical sce- 
narios in which both high magnetic fields and strong gravitational 
fields are involved, such as the magneto-rotational collapse of stel- 
lar cores, the collapsar model of gamma-ray bursts, and the evo- 
lution of magnetized neutron stars. The code is based on high- 
resolution shock-capturing schemes for solving the GRMHD equa- 
tions, which are cast in a first-order, flux-conservative hyperbolic 
form, supplemented by the flux constraint transport method to en- 
sure the solenoidal condition of the magnetic field. In addition, 
the code can handle several EOS, from simple analytical expres- 
sions to microphysic al tabulated ones. The rob ustness of the code 
was demonstrated in ICerda-Duran et all ( 1200 Sl) using a number of 
stringent tests, such as relativistic shocks, highly magnetized flu- 
ids, equilibrium configurations of magnetized neutron stars, and 
the magneto-rotational core collapse of a realistic progenitor. Al- 
though we use a nonlinear numerical code, we restrict attention to 
small perturbations and focus only on the linear behaviour of the 
oscillations. 

The fundamental variable in our simulations is the velocity 
component v"^ and the corresponding displacement Y, which we 
define through the relation 

y = Qu"^, (49) 

where a dot denotes a time derivative. Following the presentation 
of our anelastic approximation in Sec. |2] we summarize here the 
particular assumptions made in our numerical implementation: 

(i) We assume a fixed spacetime metric (Cowling approxima- 
tion). Since we consider nonrotating, stationary configurations, the 
shift vector vanishes, — 0, for the gauge choice made in the 
spatial conformal flatness approximation. 

(ii) The evolution of the continuity equation is not performed, 
since this is the main assumption of the anelastic approximation. 

(iii) Since we restrict our study to small- amplitude oscillations, 
the i;,.and ve components of the 3-velocity couple only weakly to 
the dominant component v^p and can therefore be neglected. 

(iv) In axisymmetry and with Vr = vg — 0, the poloidal part of 
the magnetic field, Br and Be, does not evolve in time and can be 
kept constant throughout the simulation. Thus, only B^ is evolved 
in time and describes torsional oscillations (axial parity). 

(v) The combination of the previous two assumptions makes it 
possible to avoid the evolution of the momentum equations for Sr 
and Se- This approach simplifies the requirement given by Eq. i fTSl l, 
as it is trivially fulfilled in our case. 

(vi) The boundary conditions at the surface of the star are as- 
sumed to be the zero-traction conditions for shear waves in a solid 
crust (in order to mimic the coupling of different magnetic field 
lines by a solid crust, see Sec. 14. 3t . For this corresponds to 
drY — 9r(Qi'*') = 0, which is equivalent to drY = if the 
initial conditions also fulfill that drY{t = 0) = 0. The boundary 
condition for B'^ is computed numerically from the condition for 
v"^ by means of the induction equation. 

With all these assumptions only the equations for S,p and B^ 
are evolved in our anelastic simulations of torsional Alfven modes. 
We choose an initial perturbation of the form 

Y{t = 0)^0 ■ y(t = O) = /(r)6(0), (50) 



where 
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Figure 3. Spatial pattern of the effective amplitude for several upper QPOs, obtained in the anelastic approximation for the reference model MNS2. Both 
antisymmetric (upper panels) and symmetric (lower panels) oscillations are shown. The thin black lines are magnetic field lines, while the thick black line is 
the stellar surface. The nearly horizontal (blue) lines are the nodal lines predicted by the semi-analytic model for standing-wave oscillations in the region of 
open field lines. 



h{e) = h2^deP2{cose) + h3^deP3{cose), (51) 

where P2 and P3 are the corresponding Legendre polynomials. The 
coefficients 62 and 63 are such that the maximum initial amplitude 
of Y throughout the star is just b2 +63. By setting 62 = the 
perturbation is symmetric with respect to the equator, and by setting 
&3 = it is antisymmetric. Any other combination has no defined 
symmetry. 

It is worth noting here that the zero traction boundary condi- 
tion implies a non-zero angular momentum flux at the surface, since 
F^{Sip) = —bipB^/W. Because of this, the star acquires a net 
angular momentum. This effect is just an artifact, arising from ne- 
glecting the presence of a magnetosphere, where part of the Alfven 



wave energy would be transmitted (the actual boundary conditions 
would then be the continuity of a nonvanishing traction across the 
surface). Notice that, antisymmetric perturbations add no net angu- 
lar momentum to the system, while symmetric perturbations intro- 
duce some amount. For small enough oscillation amplitudes (such 
as those considered here) the net angular momentum added has a 
negligible effect on the evolution of torsional modes. 

Indeed, the radial function /(r) is chosen to satisfy the bound- 
ary conditions at the surface of the star and to minimize the net an- 
gular momentum added by the perturbation. We have checked that 
the net angular momentum added by this particular choice is much 
smaller than the angular momentum introduced by the boundary 
conditions during the evolution. Furthermore, we note that due 
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Figure 4. Spatial pattern of the lower QPOs (reference model). Color coded is the oscillation amplitude, y{r, 6; /), for the fundamental lower QPO (Lq) and 
the first two overtones (Li and L2). The thin black lines show the magnetic field line structure inside the star and the thick black line represents the stellar 
surface. 

to particular choice i50\ for the initial perturbation and since the 
equilibrium magnetic field is poloidal, the initial value for B'^ is 

B-^{t = 0) = 0. 



5.1 Simulations for the reference model 

We begin our simulations by studying the evolution of the reference 
model MNS2. To this equilibrium model we add an initial velocity 
perturbation, consistent with Eq. i50l . with 02 = 03 = 5 x 10^''. 
This perturbation has no defined symmetry and its amplitude is 
small enough to be considered at the linear level. By adding this 
perturbation torsional oscillations are induced in the star which can 
be observed as oscillations of v'^ and in time. 

We perform simulations using, initially, two different numeri- 
cal resolutions, rir x ng = 80 x 40 and 80 x 80, where Ur is the 
number of cells in the radial direction and ng in the angular direc- 
tion, for a full-grid setup (9 up to 180°). We integrate in time for 
800 ms, which corresponds to about 40 times the Alfven crossing 
time, which is the dynamical timescale in our anelastic simulations. 
As expected, the amplitude of oscillations is damped throughout 
most regions in the star, due to the existence of a continuum of fre- 
quencies. At the same time, we observe persistent, long-lived QPOs 
near the pole and within the region of closed field lines, as expected 
by ou r semi-analytic model and in agreement with ISotani et al.l 
1 I2OO8I) . 

Since we do no expect global collective oscillations (normal 
modes) to be present, the analysis of the oscillations has to be 
performed locally. We denote by A{r, 9; f) the amplitude of the 
Fourier transform of y(r, 9;t) at each point in the star, correspond- 
ing to an effective oscillation amplitude at a given frequency, for a 
given, finite, simulation time. At the frequencies predicted by our 
semi-analytic model we find strong QPOs near the magnetic axis 




Figure 5. Time evolution of the maximum perturbation amplitude, 
imax(xi (color plotted, logarithmic scale), for each magnetic field line, 
X, of the reference model with a resolution of 80 X 80. The color scale 
ranges from zero (white) to its maximum value (red-black). 
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Figure 6. QPO continuum in the anelastic approximation. The logarithm of the effective amplitude, < A{x\ f) >, averaged along each individual magnetic 
field line, is shown in a frequency vs. magnetic field line plot. The central and right panels correspond to antisymmetric and symmetric initial perturbations, 
respectively, while the left panel corresponds to a mixed initial perturbation with no symmetry. The (color) scale is the same in the three panels, ranging from 
a minimum value (white) to the maximum value (black). Solid lines show the continuum of standing-wave frequencies (repeated at different overtones), as 
obtained by our semi-analytic model, for open (red) and closed (blue) magnetic field lines inside the star. 



and within the region of closed magnetic field lines. For overtones, 
the number of nodes agrees also with the predictions of the semi- 
analytic model. Once we identify the QPOs, guided by the semi- 
analytic model, we refine the determination of the QPO frequencies 
by considering the maximum effective amplitude near the magnetic 
axis and within the region of closed magnetic field lines. The val- 
ues of different extracted QPO frequencies (up to a certain over- 
tone) are displayed in Table (2] along with their relative difference 
when compared to the predictions of the semi-analytic model. Al- 
though our anelastic approach and the semi-analytic model employ 
very different approximations, the agreement in the QPO frequen- 
cies is remarkable. The frequency of the fundamental QPO agrees 
to within 4% between the two approaches, while the agreement in 
all other frequencies is at the 1-2% level. 

Fig. [3] displays the spatial pattern of the effective amplitude 
A for symmetric and antisymmetric "upper" QPOs, up to the third 
overtone. Notice that A has been rescaled by r sin 9 in order to 
correspond to a velocity component in a unit basis and not in the 
coordinate basis - in the latter the amplitude is maximum exactly at 
the pole. For all overtones, the maximum amplitude of the QPOs 
appears at the surface, at a location near the magnetic axis. As 
expected, the u!c^ QPOs are antisymmetric w.r.t. the equatorial 
plane, while the f/^^^ QPOs are symmetric. The number of nodes 
along the magnetic axis and their precise location agrees with the 
prediction of the semi-analytic model (the nearly horizontal lines 
in this Figure are the predicted nodal lines). 

Fig- El displays the spatial pattern of the effective amplitude 
A (again rescaled by r sin 9) for the "lower" QPOs, up to the sec- 
ond overtone. The amplitude of the "lower" QPOs is appreciable 
only within the region of closed field lines (in the third panel of 
Fig. [4] a nonzero pattern near the magnetic axis is a contamination 



of the finite FFT by other QPOs). As predicted by our semi-analytic 
model, the fundamental QPO, Lo, has one nodal line crossing the 
region of closed field lines, while each consecutive overtone has 
an additional nodal line, each time added in a symmetric fashion. 
Since the initial perturbation used in this simulation had no partic- 
ular symmetry w.r.t. the equatorial plane, the precise location of the 
nodal lines inside the region of closed field lines cannot be easily 
be predicted by the semi-analytic model. The resulting patterns are 
nevertheless very close to those obtained for simulations that excite 
only QPOs of a particular symmetry, i.e. to Lq^\ L[ ' and Lj^'. 

Using the magnetic field structure of the unperturbed star we 
can map the time-derivative of the displacement, Y{r, 9; t), from 
polar coordinates to the magnetic-field-line-adapted coordinates, 
y X', t)- This allows us to define the absolute value of the max- 
imum of Y along each magnetic field line, at any given time, i.e. 
|ymax(xi 01- The time-evolution of this quantity, for different mag- 
netic field lines is shown in Fig.|5] for the 80 x 80 reference model. 
Here, one can clearly observe how the initial perturbation evolves 
into long-lived QPOs near x ~ and x ~ 0.9. In contrast, the 
oscillations at intermediate points are damped in time (see Sec. 15.41 
for a discussion on the damping rates). 

In a similar way, one can map the effective amplitude 
A{r, 9; f) to .4(^, X: f) and then compute its average value along 
each magnetic field line, as 

<A{x;f)>^ [ d^A(i,x;.f)- (52) 
Jo 

The left panel of Fig.|6]shows the distribution of < A{x', /) > for 
different magnetic field lines and frequencies. One can see that the 
effective amplitude displays local maxima at the predicted QPOs. 
In addition, a broader agreement with the continuum frequency 



Alfven QPOs in magnetars 1 1 



parts for standing waves, predicted by the semi-analytic approach 
is also evident (red and blue lines overplotted). Since no particular 
symmetry was present in the initial perturbation for this simula- 
tion, both symmetric and antisymmetric upper QPOs appear. In the 
middle and right panels we show similar plots, but for antisym- 
metric or symmetric initial data, respectively. This time, only the 
corresponding sub-families of upper QPOs appear. Moreover, the 
lower QPOs (near x ^ 0.9) appear at all integer multiples, since 
the frequencies of symmetric and anti-symmetric QPOs of same 
order coincide. Fig.[6]demonstrates that (at least for low amplitude 
perturbations) the QPOs excited by symmetric and antisymmetric 
perturbations completely decouple. This may be relevant to the ap- 
pearance of only odd-integer multiples in the observed QPO fre- 
quencies in SGR 1806-20. 

5.2 Validity of the semi-analytic approach 

The semi-analytic model of Section |4] predicts that both the up- 
per QPOs and lower QPO overtones should appear at exact inte- 
ger multiples of their corresponding fundamental frequency. The 
results presented in Table |2] agree with this to within 2%, which 
is comparable to the frequency resolution of the FFT, which was 
1.25 Hz. The only exception is the fundamental upper QPO fre- 
quency, which shows a larger difference w.r.t. the expected value, 
based on the obtained values for the frequencies of the overtones, 
that cannot be explained by the finite frequency resolution of the 
FFT. We have checked that this is indeed not a numerical arti- 
fact, since the numerical frequency does not converge to the semi- 
analytic result with increasing resolution. As a result of this de- 
parture of the fundamental upper QPO frequency from the predic- 
tion of exact integer multiples one can observe an interference of 
the first antisymmetric overtone, U[~\ with three times the fre- 
quency of the fundamental QPO, (/q ^, which produces a beat 
frequency of /beat = (3f/o"' - U[^^)/2. Note that the fre- 
quency of the beat does not change appreciably with resolution. 
Measuring the beat frequency, the difference (3 Uq~^ — U\~^) can 
be computed with higher accuracy than from the FFT itself. Us- 
ing the last three "beats" in the simulation we obtain a value of 
(3 Uq~'' — (7} ') = 6.6 Hz with a frequency resolution of 0.4 Hz. 
The origin of this frequency difference with respect to the semi- 
analytic model is the short-wavelength assumption in the latter. 
However, for the fundamental upper QPO the wavelength of the 
corresponding standing-wave is twice the size of the star itself, and 
therefore our approximation becomes inaccurate. For all other over- 
tones, however, the short-wavelength approximation is valid with 
good accuracy. 

5.3 Dependence on the symmetry of the initial perturbations 

In order to check the dependence of our numerical simulations on 
the symmetry of the initial perturbation we have performed two 
more simulations of model MNS2 with a symmetric (a2 = and 
a-i = 10~^) and an antisymmetric perturbation (02 = 10^^ and 
as = 0). Since the background model is symmetric with respect 
to the equatorial plane and the perturbations are either symmetric 
or antisymmetric, we can perform the simulations using only one 
quadrant (90°), and imposing the appropriate symmetries in the 
equatorial plane depending on the initial perturbation. Therefore, 
the boundary conditions for the perturbed quantities, and B"^, 
in the equatorial plane, depend on the symmetry of the perturbation 
with respect to the equator: for symmetric perturbations, dgv^ = 
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Figure 7. Time evolution of the maximum absolute value of v^p during 
oscillations, at three representative points inside the star: at the points where 
the upper QPO (upper panel ) and lower QPO (middle panel) show their 
respective maximum effective amplitudes and at r = 0.65re and d = ■n/i 
(lower panel), a point on an open field line, far from both regions where the 
upper and lower QPOs are prominent. Three resolutions were used:, 80 X 20 
(green dots), 80 X 40 (blue dots) and 160 X 40 (red dots). Straight lines are 
fits to an exponential law for t > 100 ms. 



and B'^ = 0, and for antisymmetric perturbations, v"^ = and 
9a -B*^ = 0. We have performed simulations using 80 x 40 and 
80 X 20 grid cells in (r x 9), which correspond to the same res- 
olution as the full 180° simulations of 80 x 80 and 80 x 40 cells 
respectively. We have also performed low resolution simulations of 
80 X 40 cells in a 180° domain, checking that there is complete 
agreement with the 90° simulations with the same effective resolu- 
tion of 80 X 20 cells. 

When using either symmetric or antisymmetric initial pertur- 
bations, the spatial pattern at each QPO frequency is similar to the 
corresponding cases in Fig. [3] and the number of nodes and their 
location match as well (see also central and right panels of Fig.|6]l. 
The location of the nodes of the lower QPOs can be predicted be- 
cause of the particular symmetry of the perturbations. Similarly, for 
the lower QPOs, the spatial pattern for each QPO frequency agrees 
in the number of nodes and their location with the corresponding 
semi-analytic model, for either symmetric or antisymmetric ini- 
tial perturbations. A 90° spatial phase difference is found between 
symmetric and antisymmetric perturbations, which is expected for 
the type of initial data we use. 

We note that the time evolution of the perturbations in all the 
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different background models we use is qualitatively similar to the 
results reported using the reference model. 



5.4 Damping rates and viscosity 

According to our analysis of section|4l long-lived QPOs (upper and 
lower) are expected near the magnetic axis and inside the region of 
closed field lines, corresponding to extrema in continuum parts of 
the spectrum. In other regions, initial oscillations should be damped 
quickly, due to phase mixing. Apart from this, our simulations are 
affected by intrinsic numerical viscosity (due to finite-differencing) 
which produces an additional damping of the oscillations every- 
where in the star. We use our reference model MNS2, with an an- 
tisymmetric initial perturbation, in order to examine the effect of 
numerical resolution on the damping rate of oscillations at differ- 
ent regions inside the star, i.e. in order to examine the effect of 
numerical viscosity. We compute the evolution of this model with 
low resolution (80 x 20), standard resolution (80 x 40) and high 
resolution (160 x 40), for a half-grid setup. Fig.|7]shows the time 
evolution of the maximum absolute value of v^p during the oscil- 
lations, at three representative points inside the star: at the points 
where the upper QPO (upper panel ) and lower QPO (middle panel) 
show their respective maximum effective amplitudes, and at a point 
on an open field line, but far from the regions where both the up- 
per and lower QPOs are prominent (lower panel), specifically, at 
(r,e») = (0.65re,7r/4). 

For the two QPO families (upper and middle panels of Fig.[7]l 
we observe that the damping is significantly reduced at high res- 
olution. Especially for the upper QPOs (upper panel), already 40 
angular cells seem to be sufficient to have practically negligible nu- 
merical damping. In contrast, the lower QPOs are always damped 
in time, even with the highest resolution we use. This could be 
due, in part, to the fact that closed magnetic field lines are much 
shorter than open field lines, and therefore much higher resolution 
is needed to get a convergent result. From these results it is difficult 
to conclude what the actual damping of the lower QPOs would be, 
as even higher resolution simulations would be necessary to further 
diminish the numerical viscosity. Finally, we observe a very strong 
damping rate for oscillations at the point which is far from the two 
regions where the upper or lower QPOs are prominent (lower panel 
of Fig.|7](, and this rate does not decrease appreciably with increas- 
ing resolution. This indicates that the damping at such points must 
be due to phase mixing. 

In some magnetohydrodynamic systems, quasi-modes can be 
found, which are no rmal modes w ith a very large imaginary part in 
their frequency (see lLeviiil ( |2007|) and references therein). We have 
investigated the relation of the QPOs found here to possible quasi- 
modes in the limit of high viscosity, by adding artificial viscosity 
to our simulations. In particular, we performed a simulation of the 
reference model MNS2, with a resolution of 80 x 20 cells (in a half- 
grid setup) and an initial antisymmetric perturbation, adding artifi- 
cial viscosity by mea ns of a Kreiss-Oliger 4 th order term in the ip 
momentum equation l lKreiss&01igerlll973h . We used a large pre- 
factor in this viscosity term equal to 0.1. In this case we indeed find 
results that are consistent with global collective oscillations with 
a large imaginary part (quasi-modes). The computed real parts of 
quasi-mode frequencies and their phase show a minimal variation 
for different magnetic field lines (x) while the real part of the fre- 
quency of quasi-modes is equal to the real part of the frequency of 
the corresponding QPOs in the non-viscous case. 



5.5 Empirical relations 

In other to check the dependence of the QPO frequencies on the 
magnetic field strength and on the mass of the neutron star we also 
compute the evolution of models for which we have varied either 
the central current density (models MNS 1 and MNS3) or the cen- 
tral density (models HMNS2 and LMNS2), with respect to the ref- 
erence model MNS2. The equilibrium properties of these models 
are shown in Table[T] The simulations are performed in a half-grid 
setup with a resolution of 80 x 20 cells. Our results are summarized 
in Table [3] Frequencies computed with the semi-analytic, short- 
wavelength approximation are labeled "sa", while those obtained 
in the anelastic approximation are labeled "num". In parentheses 
we show the relative difference with respect to the corresponding 
semi-analytic value. As for the reference model MNS2, we observe 
a relative difference of several percent for the frequency of the fun- 
damental upper QPO. As we explained in Sec. 15.21 this is because 
the wavelength of the fundamental QPO is too large for the short- 
wavelength approximation to be more accurate. However, since fre- 
quencies of all overtones show a relative difference of order 1-2% 
only, with respect to the corresponding frequencies of the semi- 
analytic model, if we assume that these are exact integer multiples 
of a fundamental frequency, then we can construct a fit that yields a 
value for the fundamental frequency compatible with this assump- 
tion. These values are labeled as "fit" in Table[3] 

When varying the strength of the magn e tic fie ld we find that, 
as expected from the results of ISotani et alj j2008t) . the QPO fre- 
quencies change linearly with B. Also, the effect of changing the 
stellar mass is consistent (within a few percent) with the expansion 
in terms of the co mpactness, (M/R), used in the empirical rela- 
tions presented by ISotani et alj ( |2008|) . This agreement allows us 
to construct empirical relations for all QPO frequencie s, taking ad- 
vantage of the {M/R) expansion in'S otani et alj ( l2008h . which was 



constructed for a set of different EOS (while here we only consider 
a polytropic EOS). We note tha t in pra ctice, the larger number of 
data points used in ISotani et alj ( 1200 Sl) results in a more accurate 
empirical relation than using the few models presented here. In or- 
der to construct new empirical relations, we use the frequencies 
labeled "fit" in Table|3] 

Our empirical relations (constructed as described above) are 
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for the family of upper QPOs. Here, n indicates the order of the 
QPO, with n = corresponding to the fundamental QPO. Notice 
that the above two relations include all obtained QPOs, irrespective 
of their symmetry with respect to the equatorial plane. For the lower 
QPOs, the symmetric and antisymmetric cases have coincident fre- 
quencies and appear at all integer multiples. For the upper QPOs, 
the antisymmetric cases correspond to n = 0, 2, 4, .. (odd integer 
multiples) while the symmetric cases correspond to n = 1, 3, 5, .. 
(even integer multiples). Such empirical relations are useful for try- 
ing to interpret observed QPO frequencies and for extracting useful 
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Table 3. Frequencies (in Hz) of the lower and upper QPOs for all models. Values computed with the 
semi-analytic standing-wave approximation are labeled "sa", numerical values of the fundamental 
frequency obtained in the anelastic approximation are labeled "num" and fitted values using the 
harmonics are labeled "fit" (see text for details). In parentheses we show the relative difference with 
respect to the corresponding semi-analytic value. 



Model 


..(-)sa 

'-'o 


..(-)num 
^0 


!7l-'*^V(2n + 1) 


r sa 
^0 


r num 
^0 


L«V(n + l) 


MNSl 


3.63 


3.50 (3.6%) 


3.62 (0.3%) 


4.17 


4.12(1.2%) 


4.23 (1.4%) 


MNS2 


36.3 


35.1 (3.3%) 


36.8 (1.4%) 


41.7 


41.3 (1.0%) 


42.16(1.1%) 


MNS3 


91.2 


87.4 (4.2%) 


89.2 (2.2%) 


104.7 


103.1 (1.5%) 


106.2(1.4%) 


HMNS2 


35.4 


33.7 (5.0%) 


35.9(1.5%) 


43.3 


43.73 (1.0%) 


43.4 (0.2%) 


LMNS2 


35.2 


33.7 (4.2%) 


35.0 (0.6%) 


39.2 


38.71 (1.2%) 


39.51 (0.8%) 


SI 


21.7 


20.0 (7.8%) 


21.96(1.2%) 


25.4 


25.95 (2.2%) 


25.9 (2.0%) 



infoimation about the characteristics of the source (see the discus- 
sion in Sec.|6l(. 

We can now apply the above empirical relations to the two 
known SGR sources that exhibit the QPO behavior in their X-ray 
tail. For this , we us e the set of equilibrium models presented in 
ISotani et al which were constructed for different tabu- 

lated EOS and different masses. Within this representative set of 
models, the compactness ratio AI/R ranges from 0.14 to 0.28. If 
we assume that the observed frequencies of 30, 92 and 150Hz in 
SGR 1806-20 are produced as QPOs near the magnetic poles (up- 
per QPOs) and are odd-integer multiples of the fundamental fre- 
quency of 30Hz, then the empirical relation ( |54t restricts the dipo- 
lar component of the magnetic field to be in the range of 5 x lO^^G 
to 1.2 X lO^^'G. Similarly, if the observed frequencies of 28, 53, 
84 and 155Hz in SGR 1900-1-14 are produced as QPOs near the 
magnetic poles (upper QPOs) and are near-integer multiples of the 
fundamental frequency of 28Hz, then the dipolar component of the 
magnetic field is restricted to be in the range of 4.7 x lO^'^G to 
1.12 x lO^'^G. 

Notice that the above values refer to the strength of the mag- 
netic field at the poles. For the particular form of the dipolar mag- 
netic field, the mean value at the surface is a factor of roughly 2/3 
smaller than the strength at the pole. Therefore, for the chosen sam- 
ple of equilibrium models and taking both SGR sources into ac- 
count, the mean surface magnetic field strength of the dipolar com- 
ponent is restricted to be in the range of 3 x lO^^G to 8 x lO^^G, 
if the fundamental QPO frequencies are identified as above. Still, 
the fundamental QPO frequencies could be any integer fraction of 
30Hz in SGR 1806-20 or of 28Hz in SGR 1900-F14. For this reason, 
the range determined above is only an upper limit to the possible 
strength of the magnetic field. Notice that using the observed period 
and spin-down rate of the magnetar in the two SGR sources, one 
can estimate the mean surface magnetic field to be about 6 x 10^* G 
for SGR 1900-^14 and 8 x IO^^'G for SGR 1806-20, which is within 
the upper limit established above. 



6 DISCUSSION 

We have presented the first two-dimensional numerical simula- 
tions of axisymmetric, torsional Alfven oscillations in magnetars 
in the anelastic approximation, in which fluid modes are sup- 
pressed. In addition to the numerical simulations, we have also 
computed Alfven oscillation frequencies along individual mag- 
netic field lines with a semi-analytic approach, employing a short- 



wavelength approximation. The continuum of standing-wave oscil- 
lations obtained using the semi-analytic approach agrees remark- 
ably well with QPOs obtained via our two-dimensional simulations 
This agreement will allow for a comprehensive study of Alfven 
QPOs for a large number of different models, without the need of 
time-consuming simulations. 

Our semi-analytic approach is not restricted to a dipolar field, 
but can be extended to other magnetic field configurations, while 
our anelastic approximation does not suffer fr om the nume r ical in - 
stabilities encountered in the linear study of ISotani et al.l ( l2008h . 
Furthermore, our anelastic approach could be extended to include 
shear waves in the crust, which may prove to be important for the 
interpretation of observed QPOs. 

Our resul t s agr ee qualitatively with those obtained 
by ISotani et al.l ilOOA who worked in the linear regime, but 
otherwise used a similar computational setup as the one employed 
here. One major difference is our finding that apart from the 
anti-symmetric upper QPOs, symmetric upper QPOs also exist, 
which completes the description of the possible QPO families. 
In addition, our implementation of the boundary condition at the 
surface, motivated by the presence of a solid crust in real magne- 
tars, results in a somewhat different spectrum than in lSotani et al.l 
([2008). Specifically, QPOs generated near the axis have a nonzero 
amplitude at the surface, which would allow these oscilla tions 
to reach the magnetosphere, while in ISotani et alj ( |2008[) the 
oscillations nearly vanished at the surface. This difference in the 
applied boundary conditions at the surface results in a differ- 
ent structure of standing-wave QPOs near the pole, hence the 
quantitat ive differenc e betw een our empirical relations and those 
given in ISotani et al.l ( |2008[) . In contrast, the QPOs produced in 
the region of closed magnetic field lines are not affected by the 
boundary condition at the surface. 

Improvements of the Alfven QPO model can be achieved 
by the inclusion of a realistic crust in the nonlinear simulations 
and of correspondingly more realistic boundary conditions at the 
crust-core interface and at the stellar surface. The crust is expected 
to oscillate predominantly at the frequencies of the turning-point 
QPOs, since the magnetized core will pull at the crust coherently 
at these fr equenc i es, wi th an amplitude that is only slowly reduced 
as 1/t^/^ iLevinI ( l2007h . In addition, the crustal torsional normal 
mode oscillations may also be triggered. Regarding the edge-QPO 
that is present in our semi-analytic model for antisymmetric stand- 
ing waves, at X ~ 0.66 (at the last open magnetic field line), even 
though we could not detect the presence of a corresponding QPO 
in our two-dimensional simulations in the anelastic approximation. 
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the situation may change in the presence of a crust and with a cor- 
responding realistic boundary condition that will couple magnetic 
field lines in this area. It is thus possible that the 1 8 and 26Hz QPOs 
in SGR 1806-20 are due to crust oscillations and/or edge QPOs, but 
additional turning-point QPOs, due to a more complicated mag- 
netic field structure than considered here, could also be expected. 
The polar-parity Alfven oscillations are also of particular interest, 
as these could be excited after a n SGR event, in a ddition to the 
torsional modes studied here (see lSotani & Kokkotas (_2009) for a 
recent study of polar Alfven oscillations in relativistic magnetar 
models). 

Other effect that has been neglected so far in previous works, 
as well as in our simulations is t he presence of a superfluid phase 
in the magnetar interior. Recentlv lAndersson et al.l I2OO8,) have sug- 
gested that the effects of including this more realistic kind of mat- 
ter should influence the frequencies by a factor of a few and could 
potentially have an impact on the computation of the QPO frequen- 
cies. If this were the case, the frequencies computed in this paper 
and the upper limits for the magnetic fields proposed could change 
by a factor of a few. 

Finally, the coupling to an exterior magnetosphere has to be 
taken into account, as well as the precise mechanism by which the 
X-ray flux in the afterglow of SGR flares is modulated by Alfven 
oscillations. We are planning to address these issues in future ex- 
tension of the work presented here. 
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